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Synopsis 


The present thesis is a study of some algorithmic aspects of /5-skeletons. The 
concept of a /5-skeleton was introduced by Kirkpatrick and Radke [KR85] as a 
way of capturing the “internal shape” of a planar set of points. Indeed, we get 
an entire spectrum of shapes by varying the parameter /5. For a fixed value of /5, 
the “internal shape” is a geometric graph, obtained by joining each pair of points 
whose /5-neighborhood is empty. 

For /5 > 1, the /5-neighborhood of a pair of points Pi,Pj is the interior of the 
intersection of two circles of radius P\pipj\/2, centered at the points (1 — p/2)pi + 
(P/2 )pj and (/5/2)p* + (1 — P/2)pj, respectively. For /5 € [0,1], it is the interior 
of the intersection of two circles of radius \pjpj\/(2p), passing through p { and pj. 
We denote this neighborhood by lunep(pi,pj). The parameter /5 characterizes the 
neighborliness of a pair of points. 

Let P = {pi , Pi , • • . ,p n } be a set of n points in the plane. The /5-skeleton of P 
is a geometric graph (P, E ) such that for every p i} pj e P, the edge p^p] e E if and 
only if the /5-lune lunep(pi,pj) is empty. 

We discuss a sweepline algorithm for constructing a /5-skeleton for any p > 
0. The time-complexity of our algorithm is in 0(n 3/2 log n) for /5 > 1 and in 
0(n?/ 2 log n), for p (E [0, 1) under the metric L v for 1 < p < co. In the metrics L\ 
and Loo, our algorithm takes 0(n 5/2 log n) time for any p > 0. 
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The concept of a /3-skeleton have natural generalizations to that of /c/3-skeletons 
and additively weighted /3-skelet.ons. We have shown that the above sweepline 
method can be extended to construct these geometric graphs too. The time com- 
plexity of computing the A;/3-skeleton is in 0(/cn 3/ " 2 log n + k 2 n ) for /3 > 1 and in 
0(n r '/ 2 log n 4- n 2 k) for (3 in the range [0, 1). The time complexity of computing the 
weighted /3-skeleton is in 0(n 3/2 log n) for /3 > 1 and in 0(n 5/2 log n) for /3 in the 
range [0, 1). 

We have obtained an optimal output-sensitive algorithm for computing a /3- 
skeleton for any /3 > 2 in metric L\ and L 0 c . The time complexity is in 0(n log n + 
A;), where k is the size of the /3-skeleton. 

We have presented efficient algorithms for computing the entire spectrum of /3- 
shapes in L 2 metric. This means computing for each pair of points the largest 
/j value, /3 max , for which the /3-neighborhood is empty. Our algorithms are in 
0(n 3 / 2 log n + K) for p rnax values not less than 1 and in 0(n 2 log n + K) time 
for values greater than or equal to 0. A point that lies in the /3-neighborhood of an 
edge is a witness to this edge. The quantity K is a count of the number of times that 
the edges of the initial graph on P are witnessed for /3 = oo. We have also presented 
algorithms for computing tin; spectra of of A;/3-shapes and additively- weighted /3- 
shapes. 

We have extend the definition of /3-skeleton into higher-dimensions. We have 
shown that the size of the /3-skeleton for ft > 2 is linear in d-dimensional space, for 
any fixed d. We have presented an efficient algorithm for computing the //-skeletons 
for /3 > 2. 
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Chapter 1 


Introduction 


Shape is only operationally defined. Thus things do not have a shape the 
way Santa Claus has a red suit - Jan Koenderink in Solid Shape [Koe90] 


In many applications in image processing, computer vision, pattern recognition, 
computational morphology, and geometric modeling et al, an object is specified in 
terms of a set of points. The shape reconstruction and subsequent identification of 
this object from the set of input points turns out to be one of the main problems 
confronting researchers working in these areas. Similarly, in spatial analysis [GB78, 
HCF77], the locational identifiers of some real spatial phenomena are represented 
by a set of points. Given the input points, the problem here is to understand and 
characterize the real phenomena. 

A set of points, however, is an ambiguous representation of an object or a 
real phenomenon and, therefore, cannot be used directly in many applications. 
On the other hand, though the whole boundary can be used to define an object 
unambiguously, such a description would take up a huge amount of storage and a 
lot of processing time. To optimize the use of resources, we therefore need efficient 
methods to extract the object shape from the point set and fast methods to perform 
various operations on the reconstructed object. 
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The complexity of reconstructing the shape of an object from a set of points also 
depends on the point extraction method. For example, the problem becomes harder 
if the input points are unorganized. In many applications, the point extraction 
method is not known and input points are unorganized. So an unified method for 
reconstructing the object shape is required to handle data from various sources. 

In some applications, the input primitives are more complex then points, through 
not unorganized. For example, a well-known problem in the field of Medical Imaging 
is to reconstruct a surface from a set of parallel cross-sectional data [BS94, GKDM96, 
AJM97], 

It is known that, in many applications geometric properties alone cannot de- 
termine the shape of an object [Tou80b]. Even the exact solution of the geometric 
problems gives only partial, approximate, or heuristic solutions to various appli- 
cation problems. Shape reconstruction using non-geometric features of points is 
beyond the scope of the thesis. 

The problem of shape recovery that we have discussed so far gives rise to the 
closely-associated question of what if anything is meant by the shape of a set of 
points. The question is important because the tools available for the shape recovery 
problem that, we have discussed so far depends crucially on our answer to this 
question. 

There is no single definitive answer. One reason for this is that there is no one 
definition that is best for all applications. So various definitions are extant in the 
literature. For example, Kirkpatrick and Radke [KR85] have defined the shape of 
a set of points as the decomposition of a point set into components that identify 
“essential” points and, among these, join “essential” neighbors. As is apparent, this 
definition depends on the interpretation of the term “essential” . 

Shape descriptors are divided into two categories: those that capture, so to say, 
the “external” shape of a set of points as opposed to those that capture what is 
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called the “internal” shape [KR85]. We explore these definitions in more details in 
the following sections. 


1.1 External Shapes 

The external shape of a set of points is obtained by identifying the “es- 
sential” extreme points of the set and, among these, joining “essential” 
neighbors [KR85]. 


Among the simplest structures that capture the external shape of a planar set 
of points are the minimum bounding box and the minimum spanning circle 
of this set. The minimum bounding box is the smallest rectangle which is parallel 
to the axies and encloses all the points. This can be computed in 0(n ) time. The 
minimum spanning circle is the smallest circle that encloses all the points. This 
can be computed in linear time [Meg83]. Instead of a minimum spanning circle 
one could consider as an alternative the minimum spanning ellipse. This can be 
computed in 0(n 2 ) time [Pos81, Pos84]. 

In the direction of increasing descriptive complexity and refinement, next to the 
above shape descriptors is the convex hull of the given point set. The convex 
hull is the smallest convex set that encloses all the points. Computing the convex 
hull of a point set is one of the most well-researched problems in computational 
geometry [PS85]. Various properties of the convex hull are used in pattern recogni- 
tion applications (see survey [Tou80b]). Depending on the application, the convex 
hull may either give too detailed a shape description or may not give enough. 
Edelsbrunner et al [EKS83, EM94] overcame this lacuna by a natural generalization 
of convex hull to the notion of a-shapes. Here a is a real- valued parameter; so we 
get a spectrum of shapes by varying a. 

The a-shape of a point set is defined in terms of generalized discs. An a-disc, 
for an arbitrary real value a, is defined as follows: 



• if a > 0, it is a closed disc of radius 1/a;; 

• if a = 0, it is a closed half-plane; and 

• if a < 0, it is the closed complement of a disc of radius — 1/a. 

For a given a, the a-hull of a point set P is the intersection of all a-discs that 
contain P. For a = 0, the a-hull is exactly the convex hull. An a-shape of a point 
set is the boundary of the a-hull, with curved edges being replaced by straight edges. 
As shown in Fig. 1, a large value of a tends to produce a crude shape of the points, 
whereas a smaller value of a gives a more detailed shape. 


lA ' ’ -V& 

Figure 1: a-shape for various valuesm of a. (Taken from [Ede87].) 

The notion of a-shapes has been further generalized to weighted a-shapes 
[Ede92, VBW94]. 

An a-shape is inadequate for a point set of variable density. For example, when 
we have two clusters of different densities. Melkemi [Mel97] proposed the notion of 
A-shapes to address this situation. The A-shape of a set of points P is defined by 
means of a set A of control points as follows. 

Let i?(A,pi) be the Voronoi region of € PU A. Let A! denote the extremities 
of the Voronoi edges that separate points in set P from those in A in the Voronoi 
diagram Vor(P U A). For q in A', let Cdisk(q ) be the closed complement of the 
Delaunay disc centered at q. 
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The Ahull of P, n(A,P) = [n qeA 'Cdisk(q)] n B(A,P), where B(A,P) = 
U peS R(A,p) 

A Delaunay edge p]p] e DT(PuA) is an Aedge if there exists a Delaunay circle 
passing through p u pj, and a point in A. The Ashape of P is the set of Aedges. 
See Fig. 2. The time required to compute an Ashape is in 0(n log n) [Mel97]. 




Figure 2: Ashape of a set of points. (Taken from [Mel97].) 


1.2 Internal Shapes 

The internal shape of a set of points is the shape exhibited by identify- 
ing the “essential” internal points of the set and, among these, joining 
“essential” neighbors [KR85]. 

The most elementary internal shape descriptor is obtained by joining all closest 
pairs of points. The closest pairs of points give an approximate density of the point 
set. This can be computed in 9{n log n) time [PS85]. Another elementary shape 
descriptor is the graph obtained by joining each point to all its nearest neighbors. 
This graph is called the Nearest Neighbor Graph (NNG). The NNG of a point set 
can be computed in 0(n log n ) time [PS85]. 

At the next level of complexity, some well-known internal shape descriptors are: 
The Sphere of Influence Graph (SIG) [Tou88], The Relative Neighborhood 
Graph (RNG) [Tou80c], The Gabriel Graph (GG) [MS84], and The Delaunay 
Triangulation (DT) [Del34]. We discuss these in detail below. 
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For each point p i: let i\ be the distance to a closest point. The circle of influence 
of pi is the circle of radius r,; . centered at p { . The Sphere-of-Influence graph is 
obtained by joining each pair of points pi and pj whose circles of influence intersect 
in more than one point. This can be computed in 0(n log n) time [Tou88]. 

Two points pj and pj are relative neighbors if and only if the following inequality 
holds. d(pi,pj ) < max{d(pi,pk),d[pj,pk)\pk € P — {PiPPj}}- Intuitively, this means 
that there is no point that is closer to both pi and p 7 than they are from each 
other. The Relative Neighborhood Graph (RNG) is a geometric graph obtained by 
joining each pair of points that are relative neighbors. The RNG can be computed 
in 0(n log n) time [Sup83, Hwa90] The RNG is a subgraph of the Delaunay trian- 
gulation (DT). Given the DT of a point set, the RNG can be computed in linear 
time [Lin94, JK87]. 

Two points pi and pj are Gabriel neighbors if and only if the circle of radius 
d(pi,Pj)/2 passing through the points pi and Pj is empty. The graph obtained by 
joining all Gabriel neighbors is called the Gabriel Graph (GG). It can be computed 
in 0(n log n) time [MS84]. The GG can be computed in linear time if the Delaunay 
triangulation of the point set is given [JT92]. 

The Delaunay Triangulation (DT) is a maximal planar descriptor of internal 
structure. It, has several interesting properties [PS85]. Two points p t and Pj are 
Delaunay neighbors if there exists an empty circle passing through them. The 
Delaunay triangulation is a graph obtained by joining all Delaunay neighbors. The 
scope of applications of this geometric graph embrace fields as diverse as biol- 
ogy, crystallography, ecology, to name a few. It can be computed in 0(n log n) 
time [PS85]. 

Just as the reshapes give rise to a continuous spectrum of external shapes by 
varying the parameter a, /2-skeletons, introduced by [KR85], give rise to a continuous 
spectrum of internal shapes by varying the real parameter (5. This can be used to 
extract the boundary of an objects [ABE97]. 
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More specifically, a //-skeleton is a geometric graph, obtained by joining pairs 
of points whose //-neighborhoods are empty. The neighborhood size depends on the 
parameter //, and are of two types. 

The circle-based //-neighborhood of a pair of points (jhpPj) is defined as follows: 
for 0 > 1, it is the union of all open circles of radius 0\pipj\/2, passing through pi 
and Pj. For 0 G [0, 1], it is the interior of the intersection of two open circles of 
radius \WPj\/ (2 //), passing through pi and pj. See Fig. 3. For 0 = 1 is exactly equal 
to the Gabriel graph. Given the DT of the point set, the circle based //-skeleton can 
be computed in linear time [KR85], for // > 1, 

[3 = 2 
[3 = 1.5 

0 = 1 
[3 = 0.5 


Figure 3: Circle-based //-neighborhood (p* , pj ) for various (3 > 0. 

The lune-based [3 - neighborhood of a pair (p h pj) of points is defined as follows: 
for 0 > 1, it is the interior of the intersection of two circles of radius P\p[p]\/2, 
centered at the points (1 — ft/2 )pi + (/3/2)pj and (0/2 )pi + (1 — 0/2 )pj, respec- 
tively. For 0 G [0, 1], the lune-based neighborhood is identical with the circle-based 
neighborhood. See Fig. 4. 
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0 = oo 

0 = 3 
0 = 2 

0 = 1 


0 = 0.75 


Figure 4: Lune-based /3-neighborhood(p, , Pj ) for various /3 > 0. 


For 0 = 1 is exactly equal to the Gabriel graph and for 0 = 2 it coincides with 
the RNG. The 0-skeleton for 1 < 0 < 2 can be computed from the DT in linear 
time [Lin94, JKY] in the metric L p , provided p = 2 or 0 = 2. 


1.2.1 Other Shape Descriptors 

Another important shape descriptor which captures internal as well as external shape 
of a point set is the 7-skeleton [Vel92b, Vel92a, Vel94]. It is defined in terms of two 
parameters c 0 and c b where c 0 ,ci € [-1, 1] and |c 0 | < \ci\. The 7-neighborhood of 
a pair of points pi and pj is defined by two circles of radius d(j>i,Pj)/ 2(1 - |c 0 |), and 
d(Pi,Pj)/ 2(1 — \c x \) passing through p* and pj such that 

• if c 0 Ci < 0 : the centers of the circles lie on the same side of the line passing 
through Pi,Pj, 
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• if cqC\ > 0 : the centers of the circles lie on opposite sides of the line passing 
through pi,pj, 

• if Ci € [—1,0] : the 7-neighborhood is the interior of the intersection of these 
circles, plus the interior of the edge connecting c 0 with c l5 

• if Ci <E [ 0 , 1 ] : the 7-neighborhood is the interior of the union of these circles, 
plus the interior of the edge connecting c 0 with Ci, 


For given points pi,pj and c 0 ,Ci satisfying the above conditions, there are two 
7-neighborhoods. Two points are 7-neighbors if and only if at least one of the 
7-neighborhood is empty. The 7-skeleton is a graph obtained by joining all the 
7-neighbors. See Fig. 5 . 



Figure 5: 7(-0.15,0.3)-skeleton of a set of points. (Taken from [Vel94].) 


1.3 Organization of Thesis 


In this thesis we study some of the algorithmic aspects of lune-based / 3 -skeletons. 
In the next chapter, we discuss a sweepline algorithm for constructing a / 3 -skeleton 
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for any /5 > 0. The time-complexity of our algorithms is in 0(n 3 / 2 log n) for /5 > 1 
and in 0(n 5 / 2 log n), for /5 G [0, 1) under the metric L p for 1 < p < oo. In the 
metrics L x and Loo, our algorithm takes 0(n 5/2 log n) time for any /5 > 1. The 
concept of a /5-skeleton have natural generalizations to that of &/5-skeletons and 
additively weighted /5-skeletons. We have shown that the above sweepline method 
can be extended to construct these geometric graphs too. In chapter 3 we present 
two output-sensitive algorithms for computing a /5-skeleton for /5 > 2, in the metries 
Li and Iqq. 

In the fourth chapter we present efficient algorithms for computing the entire 
spectrum of /5-shapes. Our algorithms are in 0(n 3 / 2 log n + K) for /5 > 1 and in 
0(n 2 log n+K ) time for /5 > 0 in L 2 metric. A point that lies in the /5-neighborhood 
of an edge is a witness to this edge. The quantity K is a count of the number of times 
that the edges of the initial graph on P are witnessed for (3 = oo. We also present 
algorithms for computing the spectra of of fc/ 5-shapes and additively-weighted /5- 
shapes. 

In chapter 5, we extend the definition of /5-skeleton to higher-dimensions. We 
show that the size of a /5-skeleton for /5 > 2 in higher dimensions is linear. We also 
present an efficient algorithm for computing a /5-skeleton for /5 > 2. 


We conclude and discuss directions for further research in the last chapter. 



Chapter 2 


Computing /3- Skeletons and Their 
Relatives 


2.1 Introduction 


In this chapter, we present a sweepline algorithm [PS85] for constructing a /3-skeleton 
for any p > 0. The time-complexity of our algorithm is in 0(n 3 / 2 log n) for P > 1 
and in 0(n 5 / 2 log n), for p € [0, 1) under the metric L v for 1 < p < oo. In the 
metrics Li and L our algorithm takes 0(n 5 / 2 log n) time for any p > 0. 

The concept of a ^-skeleton have natural generalizations to that of fc/?-skeletons 
and additively weighted /3-skeletons. We show that the above method can be 
extended to construct these geometric graphs too. 

The above algorithms are based on an efficient solution to the following n-circle 
problem: Given a set C of n circles and a set P of n points, report all the circles 
in C that are empty. 

This chapter is organized as follows. In the next section, we present an efficient 


11 



12 


algorithm for the ??.-circle problem. In the third section we show how this can 
be used to compute a /^-skeleton for any fj > 0. In the following sections, we 
define the notions of a ft/5-skeleton and additively weighted /5-skeleton. Algorithms 
for computing these are presented in the fourth and fifth sections respectively; we 
conclude in the sixth section. 


2.2 The n-Circle Problem 

2.2.1 Preliminaries 

The n-circle problem we stated in the introduction can be solved as follows: Compute 
the Voronoi diagram of the query points. Let Vi denote the Voronoi region of a point 
p For each circle Cj, we find the Voronoi region Vi that contains the center of Cj. 
The circle Cj is non-empty, if the distance between pi and its center is less than its 
radius, that is pi € Cj. The time complexity of this algorithm is in 0(n log n). But 
the problem with this method is that checking the emptiness of limes is no longer 
as simple. 

Another method of solving the n-circle problem is by repeated point location 
query in an arrangement of n circles. We could use Bentley and Maurer’s algorithm 
for this purpose [BM79], which preprocesses the circles in 0(n 3 ) time for point 
location. It takes a query point, and outputs the circles that contains this point in 
0 ( log n + K) time, where K is size of the output. 

For the problem at hand, however, we know all the query points off-line. Can 
we make use of this fact to solve our problem more efficiently? The following 
observations show that the answer is yes. 

Observation 2.1 An arrangement of n circles divides the plane into a set R of at 
most ()(n 2 ) regions. 
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Proof: Let n - 1 circles divide the plane into T(n - 1) regions. The n th circle 

intersects at most n— 1 circles. These intersection points divide the boundary of 
the n th circle into at most 2(n - 1) disjoint arcs, since there are at most 2(n - 1) 
intersection points and two consecutive intersection points define an arc. Each arc 
of the n th circle divides a region into two. So the extra regions created by the n th 
circle is at most 2(n — 1). Therefore T(n) = T(n - 1) + 0(2n - 2). This implies 
that T(n) = 0(n 2 ), since T(l)= 1. | 


Observation 2.2 If a circle is not empty, then it is sufficient to find one witness 
(that is, a point lying inside the circle) for it. 


Observation 2.3 If a point is in the region r G R, then all the circles that contain 
this region are non-empty and can be pruned. 


2.2.2 An Algorithm 

The main idea of our algorithm is that instead of pre-computing the arrangement and 
then doing repeated point location, we can achieve the same effect more efficiently 
by sweeping the circles and the points simultaneously. For details on the sweepline 
paradigm we refer the reader to [PS85]. 

We shall use the following notations and definitions throughout this chapter. 
For each position of the sweepline, we define a partial order relation < a in the set 
of regions R as follows. For two regions r { and rj in R, r { < a rj if the intersection of 
the region r { with vertical line x = a lies below the intersection of the region r# with 
that line. This ordering information is maintained in the sweepline status C. An 
event point is either an intersection point of two circles, or the leftmost/rightmost 
point of a circle, or an input point. The events are maintained in a priority queue Q, 
which is initialized to contain the input points and the left-end points of the circles. 
A circle in C is said to be active if it is either empty or a witness is yet to be found. 
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We denote the set of active circles by A c . Initially, A c = C. For each region r G £, 
let C r denote the set of active circles that contain this region. 

To describe the algorithm it is enough to specify the actions at each type of 
event point as the sweep proceeds. 

Let the event e correspond to the left end point of a circle c. We insert the 
right-end point of the circle in Q and locate the region r e £ that contains this 
event point. The boundary of c splits r into three regions r x , r 2 , and r 3 , as shown 
in Fig. 6. 

Comment : It may be noted in passing that r x and r 3 are the same region r as can 
also be seen from the Fig. 6; however , we distinguish between them because in the 
sweepline status a region is identified and maintained by the two circular arcs that 
bound it from above and below. Thus r appears in two incarnations in the sweepline 
status in this case. i 

Therefore, we replace the region r by the regions r x , r 2 , and r 3 in £. Simultane- 
ously, we create the active circle lists C n , C T2 , and C r3 ; the first and third are really 
copies of C r , while the second is obtained from C T by adding to it the circle just 
encountered. If c intersects the circular arc that bounds r x from above we insert this 
intersection point in Q. Similarly, for the circular arc that bounds r 3 from below. 
These intersection points are q x and q 2 in Fig. 6. 

Let the event e corresponds to the right end point of a circle c. If this point is 
on the boundary of the regions r x , r 2 , and r 3 , then the region r 2 is deleted from £ 
and the regions r x and r 3 are coalesced into a new region r i3 (see Fig. 7). The active 
circle list C ri3 of this new region is the same as that of r x or r 3 . We insert into 
the event queue the intersection point, if any, of the bounding arcs of the region r 13 
(point q in Fig. 7). This intersection point can very well be the right-end point of a 
circle as shown in Fig. 8. 

Let the event e correspond to an intersection point of two circles c x and c 2 . If e 
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is on the boundary of the regions r l5 r 2 , and r 3 , then these regions transform into 
new regions r[, r 2 and r' z respectively. The region r[ has a new lower boundary 
vis-a-vis r\ and the same active circle list; r 2 has new upper and lower boundaries 
as compared to r 2 and an active circle list obtained from C r2 by deleting/adding 
the intersecting circles C\ and c 2 ; r [ , has a new upper boundary compared with r 3 
and the same active circle list. If the bounding arcs of r[ intersect, we insert their 
intersection point in Q. We do the same for r' z . In Fig. 9 these intersection points 
are q\ and g 2 . Notice that the latter is the right end point of one of the intersecting 
circles. 



Figure 9: Update of £ at an intersection point of two circles. 


Let the event e correspond to an input point , we locate the region r 6 £ that 
contains this point and remove all the circles in C r (see Fig. 10). The removal of 
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each such circle involves deleting it from the set of active circles A c , deleting the 
intersection points of circle c with other circles, which are present in Q, and updating 
£. 



Figure 10: Update of C at an input point. 


2.2.3 Analysis of the Algorithm 

We claim the following result. 

Lemma 2.1 The n-circle problem can be solved in 0(n 2 log n ) time, using 0(n 2 ) 

space. 

Proof: At any abscissa, the sweepline has at most 2 n intersections with the 

circles; hence at most 2n + 1 segments. Since each segment corresponds to a region, 
the sweepline status can have at most 0{n) regions. 

Any data structure used for implementing the sweepline status £ should support 
efficient implementations of the following operations: addition and deletion of a 
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legion, point location, and merging of two adjacent regions. A height-balanced tree 
meets these needs eminently. The regions are stored at the leaves and the internal 
nodes correspond to the circular arcs' that make-up the region boundaries. The 
active circle list C r of each region r in the sweepline status is also maintained as a 
height-balanced tree. 

We maintain a Boolean array Prune\\ to keep track of the status (empty /non- 
empty) of all the circles. Initially, Prune[i } is false for each i in 1 < i < n. During 
the sw r eep, if circle a is found to be non-empty we set Prune[i] to true. 

The event queue, as well as the array Prune[n ], can be initialized in O(n) time. 

We analyze separately the complexity of processing each of the four types of 
events in Q. 

Let e be an event that is the left end point of a circle c. The cost of locating 
this point in £ is in 0(log n). If r be the region that contains this event point, 
we replace it by three new regions r x , r 2 , and r 3 in £ and if the bounding arcs of 
the region r x intersect, insert their intersection point in Q. We do the same for the 
region r 3 . The complexity of this is in 0(log n) again. Creating the active circle lists 
corresponding to these regions is in 0(n), since we have to make three copies of C r 
and insert the newly encountered circle into one of these. The cost of processing the 
left-end point of a circle is in 0(n ) time. The total cost of processing the left-end 
points of all the circles is in 0(n 2 ) time. 

A similar analysis shows that the cost of processing the right-end point of a 
circle is in 0(log n). The total cost of processing the right-end points of all the 
circles is in 0(n log n). 

Let the event e be an intersection point. The cost of locating this point in £ and 
hence the three consecutive regions, r u r 2 , and r 3 , that have e is on their boundary is 
in 0(log n). The cost of replacing these regions in £ by the corresponding regions r [ , 
r' 2 and r' z is also in 0(log n). The cost of creating the active circle lists corresponding 
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to these regions is in 0(log n), since we only need to make two deletions/insertions 
into the list corresponding to r 2 to get the active circle list of r' 2 . The active circle 
lists of r[ and r 3 are same as that of r x and r 3 respectively. We insert in Q, the 
intersection points of the arcs that bound the region r[ as discussed earlier. We do 
the same for r 3 . The complexity of this is in O(log n) again. So the cost of processing 
one intersection point is in 0(log n) time. The cost, therefore, of processing 0(n 2 ) 
intersection points is in 0(n 2 log n). 

Let the event e be an input point. The complexity of locating this point in the 
sweepline status is in 0(log n). If r be the region that includes this point, we delete 
all the circles from its active list C T and set the appropriate array entries in Prune [] 
to true. 

If the point e is on the boundary of the circle c and one of the bounding arc of 
the region r is this circle c. Let the regions r x and r 2 are respectively above and 
below the point e in the sweepline status C (see Fig. 11). The circle lists of these 
regions are differ by the circle c. That is, C n — {c} = C r2 or C r2 — {c} = C ri . The 
point e is not a witness for c, but it is a witness of all other circles of the list C Tl 
and C T2 . Therefore, we prune all the circles in C n — {c}. 



Figure 11: Update of C at an input point, if it is on the boundary of a region. 
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The cost of deleting a circle c* is in 0(n log n). We argue as follows. As 
observed earlier, this involves removing from the event queue its intersections with 
other ciicles that are present in Q as also its right end-point; it also involves removing 
this circle from the active list of all the regions in £ in which it occurs. The first can 
be done in constant time by the following artifice. We merely set Prune[i] to true. 
While processing an intersection point q or a right end point we have to take the 
additional care of checking suitable array eirtries in the array Prune]] to decide if 
q needs to be processed. The cost of the second is in 0{n log n ). The adjustments 
required to £ is in 0(log n) since the removal of its upper /lower arc causes the 
regions above and below this arc to be merged into one. In the course of the sweep 
at most n circles are pruned, and thus the amortized cost is in 0(n 2 log n). 

Hence the n-circle problem can be solved in 0(n 2 log n) time using 0(n 2 ) space. 

I 


The complexity of the above algorithm is worse than the brute-force algorithm 
which tests a point for inclusion in each of the n circles! For the complexity of this 
is in 0(n), and therefore in 0(n 2 ) for n points. However, if we look at the above 
algorithm very carefully we notice that the problem lies in considering all the circles 
in one go. Tin* processing time for the n circles is high and dominates the algorithm. 
The following lemma shows how to balance the processing time for circles and points 
to improve the time complexity. 


Lemma 2.2 Given a set C of n circles and a set P of n points, all the empty circles 
can be computed in 0(n 3 / 2 log n) time using 0(n) space. 


Proof: We divide the problem into n/m subproblems, each subproblem consisting 

of to circles and all n input' points. For each of these subproblems the above 
algorithm takes 0(m 2 log m + n log rn) time. Therefore, the total time required 
to solve all the n/m problems is 0(nm log to + n 2 /m log to) time. The minimum 
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time of 0(n 3 / 2 log n) is achieved for m = y/n. The space required is in 0(n), since 
each subproblem consists of y/n circles. Hence the claim of the lemma. ■ 


2.3 /3-Skeleton 


In this section, we generalize the above algorithm to compute the /5-skeleton for 
/?>!• 

It is well-known that the Delaunay triangulation (DT) [PS85], is a super-graph 
of a /5-skeleton for j5 > 1. Therefore, to compute the /5-skeleton, for some /5 > 1, we 
construct DT(P) and prune the edges whose /5-lunes are not empty. The Delaunay 
triangulation can be constructed by any standard algorithm in 0(n log n) time. 
What remains is to find a way to prune the non-empty lunes. 

It is easy to extend the solution to the n-circle problem, if circles are replaced by 
lunes. The same observations have been made by Su and Chang in [SC91b]. In fact, 
we can consider any set of convex objects with constant-size description, any two of 
which have a bounded number of intersections. With lunes, the event structure is 
slightly more complicated. Unlike in the case of circles, each lune can have up to 
four events associated with it. Two of these correspond to positions where a lune 
joins the sweepline and leaves it. The other two correspond to the intersections of 
the circles that define a lune (i x and 12 in Fig. 12). 

Due to this event structure the region in the sweepline status that reflects its 
membership in L can be defined by three different pairs of arcs. For example in 
Fig. 12 these are the pair {a, 6}, {a,c} and {c,d} consecutively. The rest of the 
details carry over and we have the following claim: 

Theorem 2.3 For p > 1, the p-skeleton of P can be constructed in 0(n 3/2 log n) 
time using 0{n) space. 
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Figure 12: Initial events of a lime. 


Similarly, the /3-skeleton of P, for /3 € [0, 1), can be constructed by starting with 
the complete graph on P as the super graph and pruning the non-empty lunes. The 
following theorem can be proved along exactly the same lines 


Theorem 2.4 The (3 -skeleton of a set P of n points in the plane, for f} G [0,1), 
can be constructed in 0{n 5/2 log n) time using 0(n 2 ) space. 


Proof: We divide 0{n 2 ) lunes corresponding to the edges in the complete graph 

into n 2 /m subproblems. Each subproblem consists of m lunes and all n input 
points. For each of these subproblems the above algorithm takes 0(m 2 log m + 
n log m) time. Therefore, the total time required to solve all the n 2 /m problems is 
0(n 2 m log m + n 3 /m log m) time. The minimum time of 0(n 5 / 2 log n) is achieved 
for m = y/n. The space required is in 0(n 2 ), since we are start with complete graph. 
Hence the claim of the lemma. I 


2.4 k/3- Skeleton 

We can generalize the notion of a /3-skeleton by relaxing the emptiness criterion and 
allow the interior of a lune to contain up to k - 1 points; we get what is called a 
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A', 5-skeleton. The geometric graphs that we get for different values of /?, 1 upwards, 
are all subgraphs of what is called the A-Delaunay graph. 

The A-Delaunay graph is a generalization of the Delaunay triangulation on n 
points. The triangulation in this case can be characterized this way: a triangle 
belongs to the A-Delaunay graph if its circumcircle contains at most A - 1 points. 
YVe can obtain it from the Ath-order Voronoi diagram as follows: an edge piPj belongs 
to the A-Delaunay graph if a segment of the bisector of p* and pj is in the Ath-order 
Voronoi diagram [SC90]. 

Algorithms for computing A, 5-skeletons for special values of (3 are extant in the 
literature. Su and Chang [SC90] has shown how to compute the A-Gabriel graph 
(,5 = 1) in 0(A 2 n log n) time. For fixed values of A, Chang et al. [CTL92b] presented 
an 0(n 2 ) algorithm to compute the A-Relative Neighborhood Graph (/ 3 = 2). Later, 
Su et al. improved this to 0(n 5 / 3 log n) [SC91b]. An interesting application of the 
A-RNG has been made to the Euclidean bottleneck matching problem , by Chang et 
al [CTL92b]. In this problem a matching has to be found among n points in the 
Euclidean plane that minimizes the longest edge. 

To compute a A/?-skeleton for > 1, we start with the A-Delaunay graph of P 
and prune those edges whose lunes contain more than A — 1 points. It turns out 
that the above n-circle algorithm can be adapted to this case also. The additional 
work involved is to maintain a counter for the number of points inside each lune. 
We increment this count for the appropriate set of lunes whenever an input point 
is encountered during the sweep. Whenever the counter value of any lune becomes 
A, we prune that lune from further processing. The extra time required to maintain 
all the counters is in 0(A 2 n), since there are 0(A(n - A)) counters (the A-Delaunay 
graph has 0(A(n - A)) edges [SC90]) and each counter is incremented at most A 
times. We have thus proved the following theorem. 

Theorem 2.5 The k(3-skeleton of a set P of n points in the plane for f3> 1 can be 
constructed in 0(An 3 / 2 log n + A 2 n) time. 
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It is quite easy to compute a k /5-skeleton, for values of f3 in the range [0,1). We 
need to start with the complete graph on P, and prune edges which are not in the 
/c /3-skeleton. The additional cost of this is reflected in the following result. 


Theorem 2.6 The k/3-skeleton of a set P of n points in the plane, for (3 <~ [0, 1) 
can be constructed in 0(n 5 / 2 log n + n 2 k) time. 


2.5 Weighted Neighborhood Graphs 


In this section, we discuss various additively weighted neighborhood graphs. We 
assign a positive weight W{ to each point p, € P. The weight u>i expresses the power 
of the point pi in the sense that each weighted point p z - can be treated as a circle 
Bi , centered at pi, with radius equal to the weight w,. Thus the weighted distance 
d w (pi,x ), between a point pi € P and an arbitrary point x, is d(pi, x) — Wi if x is 
outside Bi , and zero otherwise. The weighted distance d w {p i7 pj), between any two 
points Pi,Pj £ P is the distance between Bi and Bj. Let D denote the set of circles 
Bi. We discuss the details of a few individual special cases 

2.5.1 Weighted Gabriel Graph (WGG) 

The weighted Gabriel neighborhood, N G (Pi,Pj), of two points p i7 pj e P that 
have disjoint circles B { and Bj is an open circle of radius d w (pi,pj)/ 2, whose center 
is on the line segment pfpf and is tangent to the circles Bi and Bj (see Fig. 13). If 
Bi and Bj intersect, then N G (pi,Pj) is empty. 

The weighted Gabriel graph of a set of points P is a geometric graph (P, E) 
such that for every Pi,Pj € P, the edge pipf € E if and only if N G (pi,pj ) is not 
intersected by any circle Bk e D — {Bi,Bj}. 
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The weighted Delaunay graph is a dual of the weighted Voronoi diagram [AE84] 
and can be computed in 0(n log n) time [Mir93]. An edge pipj belongs to this 
graph, if there is a common boundary between the Voronoi regions of pi and pj. 
The following result and its proof are from Mirzaian [Mir93]. 


Lemma 2.7 The weighted Gabriel graph is a subgraph of the weighted Delaunay 
graph if the set of circles D is pairwise disjoint. 


Proof: Let the edge pip] € WGG. The Gabriel neighborhood of pi and pj, 

Na(pi,Pj), is not intersected by any circle in D. Let the point C PiPj be the center of 
the Gabriel Neighborhood Nc(j?i,Pj). The circles Bi and Bj are at an equal distance 
from C PiPj . Moreover, no circle B k D — {Bi,Bj} is closer to C PiPj than the circles 
Bi and Bj. From the definition of a Voronoi diagram, the point CpiPj is on the 
bisector of the sites Bi and Bj. Thus by definition pjjpj e WDG [Mir93]. | 


2.5.2 Weighted Relative Neighborhood Graph (WRNG) 

The weighted relative neighborhood, NRNG(Pi>Pj)> of tw0 points p^pj € P that 
have disjoint circles B { and Bj is the interior of the intersection of two circles, of 
radii d w {jp u pj) + w { and d w (pi,Pj) + w jt centered respectively at the points p { and 
Pj (see Fig. 14). If B { and Bj intersect then N RNG {puPj) is empty. 
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The weighted relative neighborhood graph of P is defined as follows. For every 
PiiPj € P, the edge pip] is in the graph if and only if the following inequality holds. 
d w (puPj) < max{d w (pi,p k ),d w (pj,p k )}, for all p k e P - {puPj} (see Fig. 14). Or 
equivalently, if the lune is not intersected by any circle in D — {B i: Bj}. 

Mirzaian [Mir93] made a very interesting application of the weighted RNG to 
the problem of Minimum Weight Euclidean Matching. In this problem we are given 
2n points in the Euclidean plane and are required to match them into n pairs such 
that the sum of the n distances between the matched pairs is minimum. 


2.5.3 Weighted /^-Skeleton 

For all pi,pj € P such that J5,- and Bj are disjoint, the weighted lune-based neigh- 
borhood, lune/ 3 (pi,pj), is a region lying between the circles Bi and Bj and is defined 
as follows: for (5 > 1, it is the interior of the intersection of two circles of radii 
pd w {pi,pj)/2 + Wi(f3 - 1) and /3d w (pi,pj)/ 2 + Wj(p - 1), whose centers are on the 
line joining pi and pj, and which touch externally the circles Bj and Bi respectively 
(see Fig. 15). For (3 e [0,1], it is the interior of the intersection of two circles of 
radius d w (pi,pj)/2/?, passing through p\ and p'- which are a pair of closest points 
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on the circles Bi and Bj (see Fig. 15). If the circles Bi and Bj intersect, then the 
lunep(puPj) is empty for all 0 > 0. 

— : — 0 = 3 

0 = 2 

0 = 1 

Vi 

0 = 0.75 


Figure 15: Neighborhood of weighted 0-skeleton. 

The weighted 0-skeleton of P in this case is a geometric graph (P, E) such that 
for every Pi,Pj € P, the edge pip] G E if and only if lunep (pi,Pj) is not intersected 
by any circle C k € D - {B h Bj}, where D is the set of circles B { . 

When 0 = 1, lune 0 (pi,Pj ) is exactly the Gabriel neighborhood N G (puPj), and 
when 0 = 2, we get the weighted relative neighborhood, N RNG (Pi,Pj)~ 

Lemma 2.8 For a pairwise disjoint set of circles D and 0 > 1, the weighted 0- 
skeleton is a subset of the WDG. 



Proof: This lemma follows from the fact that for 0 = 1, the weighted 0-skeleton 

is exactly the weighted Gabriel graph, which is a subgraph of WDG from Lemma 2.7. 
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2.5.4 Algorithm 

In this section, we discuss an algorithm for computing a weighted /3-skeleton for 
/3 > 1. We assume that the circles in D are pairwise disjoint. 

As we mentioned in the previous section, the WDG is a super-graph of a W(3- 
skeleton for /3 > 1. The \\ DG can be constructed in 0(n log n) time [Mir93]. 

So we just have to prune the lunes that intersect any circle B{ £ D. For the 
sake of simplicity, we consider circles instead of lunes. As we observed earlier, the 
generalization to lunes is pretty straightforward. The problem can be stated with a 
bit of color as follows: 


Given a set of O(n) blue circles C and a set of n red circles D, find all 
the blue circles that are not intersected by any red circle. 


We follow closely the basic circle-point algorithm to design an algorithm for 
this case. Most of the details carry over, except when in the course of the sweep 
we meet a red circle we locate its left-end point in a region of the sweepline status 
and check for intersection with the circles which correspond to the upper and lower 
arcs that bound this region. If a bounding arc is intersected, the corresponding 
circle is removed from the active list of blue circles, and the arc is deleted from 
the sweepline status. This latter action causes two adjacent region to be merged 
into a new one and we repeat the above sequence of steps for this region. The 
process terminates when we get a region whose bounding arcs have no intersection 
with the red circle. Now we insert the red circle into the sweepline status because 
there can be yet more intersections with blue circles that have not been swept yet. 
However, such an intersection can occur only after a blue-blue intersection has been 
encountered during the sweep (see Fig. 16). So, whenever a red-blue intersection 
is encountered, prune the blue-circle corresponding to this intersection. A routine 
complexity analysis now gives the following result. 
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Figure 16: Intersection between the red circle Bk and the blue circle Cj is detected after blue-blue 
intersection tij. 

Theorem 2.9 The weighted (3-skeleton of a set of n points in the plane for (3 > 1 
can be constructed in log n) time. 

Theorem 2.10 The weighted (3- skeleton of a set of n points in the plane for p. e 
[0, 1) can be constructed in 0(n 5 / 2 log n) time. 

2.6 Conclusion 

In this chapter, we have presented an efficient sweepline algorithm for the n-circle 
problem, and used this to compute /5-skeleton, fc/5-skeleton, and weighted /5-skeleton 
under the metric L p for 1 < p < oo. In the metrics L\ and L the size of 
these geometric graphs are in 0(n 2 ) for any p. We can compute these graphs in 
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0(n 5 / 2 log to) time by starting with the complete graph and pruning the non-empty 
lunes, following the technique of this chapter. 

The algorithms presented in this chapter are not optimal. It would be intersect- 
ing to find out more efficient algorithms for computing these geometric graphs. 



Chapter 3 


An Output- Sensitive Algorithm 
for Computing the /3-Skeleton 

3.1 Introduction 


In the previous chapter we have presented an efficient algorithm for computing the 
/3-skeleton in any metric L P) for 1 < p < oo. In the metrics L x and L^, a /3- 
neighborhood is a simple rectangle compared to the complex lune in metric L p for 
1 < p < oc. Now the question is can we make use of this property to design more 
efficient algorithms? The answer to this question is yes. We present two output 
sensitive algorithms for computing a /3-skeleton in the metrics L x and L for any 
/3 > 2. The algorithms are in 0(n 3 / 2 log n + k) and 0(n log n + k) respectively, 
where k is size of the /3-skeleton. For the rest of this discussion we shall confine 
ourselves to the L <*, metric because rotating by 45° and scaling by a factor of \/2, 
the L x metric is transformed into L <». 
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3.2 Size of a /3-skeleton 


In the metric L p for 1 < p < oo, a /3-skeleton is a subgraph of the Delaunay 
triangulation [Tou80c] for /3 > 1. This implies that the size of the /3-skeleton is in 
0(n), for /3 > 1. But in the metric , the size of a /3-skeleton is in 0(n 2 ). We 
show this by constructing an example. 

Let us place n points on each vertical edge of a square as shown in Fig. 17. The /3- 
neighborhood of each edge pq , where p £ {p h , p h , . . . , p in } and q £ {p h ,p h ,..., p jn } , 
is empty for any /3 > 1. Therefore, the size of a /3-skeleton for any /3 > 1 is in 0(n 2 ). 



3.3 Algorithm 


The /3-skeleton can be computed in 0(n 3 ) time by starting with the complete graph 
as a super graph and pruning the non-empty Junes by a brute force method. This can 
be improved to 0(n 5 / 2 log n) time by the sweepline method described in previous 

chapter. 

In the metrics L\ and L oo, the /3-lune is a rectangle (see Fig. 18). Therefore, the 
emptiness query of a lune can be decided efficiently in 0(log ti) time using range 
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tree data structure [PS85]. Hence a /3-skeleton can be computed in 0{n 2 log n) 
time. By taking a closer look into the geometry of the problem, we can improve on 
this complexity further. 


P = 2 




/3 = 3 



Figure 18: The ^-neighborhood in metric L for /3 > 2. 


The main idea is to find a smaller super graph than the complete graph on P. 
We adopt Yao’s [Yao82] region approach for this. 









35 


Consider a point pi £ P. We divide the space around p± into eight regions [Yao82] , 
by drawing four lines through p % that make an angle of 0°, 45°, 90°, and 135° respec- 
tively with the x-axis. These regions are open on one side and closed on the other. 
For example, R x is closed on the a;-axis and open on the diagonal, and R 2 is closed 
on the diagonal and open on the y-axis. The other regions are similarly defined. 
These regions are called narrow regions , because for any two points Pj,Pk in a region 
dist(pj,p k ) < max{dist(pi,pj), dist(p i: p k )}. Let the m th region, for 1 < m < 8, be 
labeled Rm(pi) as shown in the Fig. 19. 



Figure 19: Narrow regions in metric L^. 

We connect each point pi to all its nearest neighbors in each region RmiPi), for 
1 < rn < 8. The graph so obtained is called the Geographical Neighborhood Graph 
(GNG). The following lemma is easily verified. 

Lemma 3.1 A (5 -skeleton for /5 > 2 is a subgraph of the GNG. 

Proof: Prom definition we know that any /5-skeleton for (3 > 2 is a subgraph of 

the /5-skeleton for (3 = 2. So it is sufficient to show that the /5-skeleton for /5 = 2 is 
a subgraph of the GNG. 

Let assume that the edge pip] € /5-skeleton for (3 — 2 . So, lune{pi,pj) is empty. 
This /5-lune is defined by two circles of radius d{p h pj), centered at p { and pj. 
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Without loss of generality we assume that pj e Rifa). Let pfp] £ GNG. This 
implies that there exist a point p k 6 Ri{pi) such that d{p u p k ) < d(pi,pj). That is, 
p k is inside the circle of radius d(p uPj ), centered at Pi . 

The definition of the narrow region d(pj y p k ) < max{d(pi,pj), d(pi,p k )}, and 

d( PuPk ) < d(pi,pj) implies d(p jt p k ) < d( Pi , Pj ). That is, p k is inside the circle of 
radius d(pi,pj), centered at Pj . 

That is, p k e lunep =2 {pi,pj ), which is a contribution. Therefore, p^pj 6 GNG. 

Hence the claim of the lemma. | 

Now the question is, for a given point how efficiently can we find its nearest 
neighbors in each region? We consider the region fyfa) for a point p { e P. The 
remaining regions R 2 , # 3 , • • • , Rs can be handled in the same way. 

We sort all the points on their x-coordinate as primary key and y-coordinate 
as secondary key. The sorted list of points is kept in an array. For each point pi, 
we find a nearest neighbor pj (in Ri( P i)) if it exists in 0(n log n) time using the 
divide and conquer approach of Guibas and Stolfi [GS83]. The rest of the nearest 
neighbors can be found as follows. 

We compute the intersection points S and N of the boundary of the region 
Ri (jPi) with the vertical line passing through pj (see Fig. 20). We locate the points 
S and N in the sorted array of points. All the points between N and S are the 
nearest neighbors of Pi in the region R\( P i), since points are sorted by y-coordinate 
as the secondary key. We denote these by p ^ , P j 2 , . . . , P j k . 


3.3.1 Pruning edges 

Our algorithm for pruning the edges whose lunes are non-empty is based on the 
following key observation. 
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Figure 20: All nearest neighbors of p, in R^. 

Observation 3.1 Let py be a nearest neighbor of Pi in the region Ri. Then for any 
P > 1, the lunep(pi,pji) is a bounded rectangle, contained within an infinite strip, 
bounded by vertical lines passing through the points pi and py . 


Let ai and bi be points (of P ) inside the strip such that lines through these, or- 
thogonal to the bounding edges of the strip, define the largest rectangle < r\, r 2 , r 3 , r 4 > 
contained in the strip whose interior is empty as shown in the Fig. 21. Note that 
it’s possible that one or both of these points do not exist. 

Clearly, all the lunes contained in this rectangle are empty. For p > 2, an easy 
decision algorithm is based on the following observation. 


Observation 3.2 The regions lunep{pi,pjf), lunep(pi,pj 2 ), . . lunep(pi,pj k ) form 
a sequence of nested neighborhoods, with lunep(pi,Pj l+1 ) C lunep{p i ,pj l ), where 1 < 
l<k- 1. 

Here is the algorithm. Starting with the smallest one, lunep(pi,p jk ), we look for 
the first non-empty lune. If none is found, all the edges between Pi and its nearest 
neighbors are in the /9-skeleton. Else if lunep(jpi,Pj m ) is the first non-empty lune 
then the edges joining pi to pj m . . -Py are pruned. 
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Obviously enough, the time required for this procedure is proportional to the 
number of non-empty lunes. 

So all that remains to be done is to show how to compute the points a { and 
Let us do this for the point b{\ the other point < 2 j can be computed similarly. Since 
we would like to do this for every point in P, the following abstract problem can be 
formulated. 


Given a set I of n horizontal intervals and a set P of n points, for each 
interval k € /, find the closest point b { 6 P, if any, that lies below it 
and within the infinite strip determined by the vertical lines that pass 
through the end points of k 


Caveat: It just happens that for our problem points and intervals are not 
distinct; each point is the left end-point of an interval. 

We outline two different solutions to this problem: the first based on the familiar 
sweepline technique, and the second based on the range-tree data structure. 


3.3.2 Sweepline Method 

We sweep the plane from top to bottom. At any instant, the sweepline status £, 
contains a relative ordering of the end-points of the intervals that have been swept 
by the sweepline and whose closest point is yet to be found. Two adjacent end- 
points in £ represent an elementary interval. With each elementary interval El, we 
associate a list of intervals C E i that contain this elementary interval. 

The event queue consists of the points of P in y-sorted order. 

Each stop in the event queue triggers off two events. The first corresponds to 
the point p. We locate the elementary interval El in £ that contains p. The point 
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p is closest to all the intervals in Cei ■ We delete these intervals from the sweepline 
status and adjust the lists of all the elementary intervals. 

The second event corresponds to the interval l whose left end-point is p. Let q 
be its right end-point. The point p divides El into two intervals Eli and EI 2 ■ We 
delete El and insert the elementary intervals Eli and EI 2 into £. We also create 
the lists Ceii and Cei 2 - We then do the same for q. For each elementary interval 
El' that is completely inside l, we update Cep by inserting l into it. 

We claim the following result. 


Lemma 3.2 The running time of the above algorithm is in 0(n 2 log n); the space 
complexity in 0(n 2 ). 


Proof: At any time, the sweepline status £ contains at most In end-points, 

corresponding to n intervals. We implement £ by a height-balanced tree that stores 
the end-points at the internal nodes and the elementary intervals at the leaves of 
the tree. 

For each elementary interval El G £ we maintain in a height-balanced tree the 
set of intervals Cei that cover El. 

The event-queue is implemented by an array. It is initialized at the beginning 
of the sweep in 0(n log n ) time. When the sweepline encounters an event, we first 
process the point p corresponding to this event. We locate the elementary interval 
El € £ that contains p. The time for this is in 0(log n). Since the point p is closest 
to all the intervals in Cei, these are deleted. If l' is one such interval, we delete 
its left-end point from £ in 0(log n) time. This causes two adjacent elementary 
intervals to be merged. The time for this is in 0(log n) . The right-end point of 
l' is similarly deleted. We also delete l' from the interval list of each El, which is 
in between the left and right end points of l' in £. This can be done in 0(n log n) 
time, since there can be at most 0{n) elementary interval in sweepline. The cost of 
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deleting an interval is thus in 0(n log n ) time. If we charge this cost to the interval 
itself, the total cost of deleting all the intervals is in 0(n 2 log n). 

We next pi o cess the associated interval l. In 0(log n ) time we locate the 
elementary interval El e £ that contains its left-end point a and insert it in C 
in 0(log n) time. This means splitting El into two intervals EI X and EI 2 , deleting 
El and inserting EI X and EI 2 into L. All this can be done in 0(log n) time. 

The interval list C Eh is the same as C E i , while C E i 2 = C EI U {/}. The cost of 
creating C E i 2 is in 0{n). 

% 

We repeat the above sequence of actions for the right end-point of l. 

Further, for each elementary interval EF that is contained in l, we update C Er 
by inserting l into it. This can be done by traversing the sweepline C from the 
elementary interval, containing left-end point to the elementary interval containing 
right-end point. The cost of all these operations is in 0(n log n) time. 

Thus the the total cost of processing n intervals is in 0(n 2 log n) time. 

Hence this problem can be solved in 0(n 2 log n) time. | 

The complexity of the above algorithm is prima facie worse than that of the 
algorithm which, for each given interval, tests by brute-force which of the n points 
is closest to this interval. The complexity of the latter is in 0(n 2 ). However, if 
we look at the above algorithm very carefully we notice that the problem lies in 
considering all the intervals in one go. The cost of maintaining the sweepline status 
of n intervals is high and dominates the algorithm. The following lemma shows how 
the complexity can be improved by balancing the processing time for intervals and 
points. 

Lemma 3.3 Given a set I of n horizontal intervals and a set P of n points, for 
each interval l{ € I , the closest point bi G P, which is within the vertical infinite 
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strip bounded by the lines passing through the end points of k and below k can be 
found in 0(n 3 / 2 log n) time using 0(n) space. 


Proof: We divide the problem into njm subproblems, each subproblem consisting 

of m intervals and all n input points. For each of these subproblems the above 
algorithm takes 0(m 2 log m + n log m ) time. Therefore, the total time required 
to solve all the n/m subproblems is in 0(nm log m + n? jm log rn). The minimum 
time of 0(n 3 / 2 log n ) is achieved for m = y/n. The space required is 0(n), since we 
have to consider yjn intervals at a time. Hence the claim of the lemma. | 


3.3.3 Range-tree Method 

The complexity can be further improved by using the range-tree data structure [PS85]. 
For the sake of completeness we briefly describe the method. We preprocess the set 
of points such that for a given a query interval k, the corresponding point b{, if it 
exists, can be located quickly. The range-tree is a rooted binary tree which allows 
us to search along two orthogonal directions, say x and y, simultaneously. It has 
two levels, the primary and the secondary. The primary level is useful for searching 
along the ^-direction, while the secondary level for searching in the y-direction. 

Each non-root node v of the range-tree T corresponds to a set points, lying 
between the abscissae L(v) and R(v). These points are stored in an array, sorted 
by their y-coordinates, in the secondary structure Y(v). The range R(v) — L(v) 
is divided into two sub-ranges such that each sub-range contains at least half the 
number of points of v. The left and right ranges correspond to the nodes LSON(v ) 
and RSON(v) respectively. The root of T corresponds to whole range, and the 
secondary level at the root node is a threaded binary tree which contains all the 
points, sorted by their y-coordinates. 
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Foi each non-leaf node v, and for each point p of y(u), there is a pointer 
LBRIDGE{p ) (respectively RBRIDGE(p )) to a point p' oiY(LSON(v )) (respec- 
tively p" of Y ( RSON (v))), which is p itself or or a point that lies immediately below 
P- 


The range-tree can be constructed in 0(n log n) time in a bottom-up manner. 
The space required is in 0(n log n) since there are 0(log n) levels and each level 
requires O(n) space. 

The application to our problem is quite simple. We narrow down the z-range in 
which a query interval l { lies. Each time we do this, we update the point p that is 
closest to the supporting line of k in the current range to the point that is closest in 
the new range. For this latter work we use the secondary structure stored at each 
node. 

For a given query interval k, the details of the search mechanism is as fol- 
lows. Find the point p which lies immediately below the line passing through l { by 
searching in Y(root) in O(log n) time. For each node v, if k is completely within 
the range of LSON(v ) (respectively RSON(v)), we search in the sub- tree rooted 
at LSON(v ) (respectively RSON(v)). Otherwise, search in both the off-springs. 
We consider searching in the left off-spring, as the right off-spring can be handled 
similarly. While proceeding to next level we can get the point which is immediately 
below the line passing through k in Y(LSON(v)), by using the pointer LBRIDGE. 
If the range of any node is completely within the given query line segment l, the 
points which are immediately below l are candidates. Of these, the required point is 
the one that is closest. See Fig. 22 and Fig. 23. The query time is thus in 0(log n) 
time. 

We summarise our results in the following lemma. 


Lemma 3.4 The algorithm for computing a (3-skeleton for P >2 in region R x takes 
0(n log n) time. 



Figure 22: partition into standard intervals. 


The above algorithm can be applied to the remaining regions by using an 
appropriate linear transformation on the input points to get the complete /3-skeleton 
for the given 0 > 2. Reporting the edges of the /3-skeleton takes time proportional 
to the size of the output. Hence follows the theorem. 

Theorem 3.5 For 0 > 2, a 0-skeleton can be computed in 0(n log n + k) time, 
where k is size of the output. 


3.4 Conclusion 

In this chapter we have described an output sensitive algorithm for computing a 
/3-skeleton in the L «, metric for 0 > 2. The complexity of the algorithm is in 
0(nlogn + k). It would be interesting to find an output-sensitive algorithm for 0 
in the range [1, 2). 








Chapter 4 


Computing the /3-Spectrum and 
Their Relatives 


4.1 Introduction 


In this chapter we propose efficient algorithms for computing the entire /3-spectrum 
of a finite planar point set in the L 2 metric. This means computing for each pair of 
points the largest /3 value, /3 max , for which the /3-neighborhood of the pair is empty. 
The motives are twofold: algorithmic and applied. From an algorithmic point of 
view this is an interesting problem for which no efficient solution was known before. 
And from an applied perspective this has rich potential for the problem of character 
recognition since the spectra of a set of points can be construed as its signature. 
Moreover, any /3-skeleton can be extracted from the /3-spectrum in linear time. 

Our algorithms are in 0(n 3 / 2 log n+K) time for /3 max > 1 and in 0(n 2 log n+K) 
time for p max > 0. A point that lies in the /3-neighborhood of an edge is a witness 
to this edge. The quantity K is a count of the number of times that the edges of 
the initial (Delaunay or complete) graph on P are witnessed for (3 = 00 . 
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The rest of the chapter is organized as follows. In the next section we present 
algorithms based on the techniques of range searching and sweepline. In the two 
following sections we discuss respectively how to compute efficiently the &/5-spectrum 
and the additively weighted /5-spectrum. We conclude in the fifth and final section. 


4.2 /3-spectrum 


In this section we propose two different algorithms for computing the /5-spectrum 
of P. 

The main idea of both the algorithms is to start with a sufficiently large super 
graph on P and compute the witness set W e for each edge e, setting /5 = oo. For 
/5 > 1 we start with Delaunay triangulation on P since it is a super-graph of the 
Gabriel graph (corresponds to /5 = 1). For /5 > 0, we start with the complete graph 
on P. We discuss only the first case as the second case can be handled similarly. 

We begin with the key observation that all the witnesses of an edge are contained 
in the lune that corresponds to 0 = oo. This an open infinite strip bounded by the 
lines normal to the edge at its end-points. We assume without loss of generality 
that no edge is horizontal, since all such edges can be handled separately in a 
preprocessing step in 0(nlog n) time. 

Let T be a transformation that maps a point p, with coordinates (a, b ), into the 
line T p with equation y = —ax -1- b and a line l with equation y = mx + c into the 
point Tj with coordinates (— m, c ) (see Fig. 24). This is known as a pont-line duality 
transform. 

Under this point-line duality transformation T, an infinite strip in the primal 
plane maps to a vertical line segment in the dual plane. So a point p E W e if and 
only if T P intersects T e (see Fig. 25). The problem of computing W e for an edge e, 
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T 

• -e > 

P 

Figure 24: Geometric dual transformation T 

therefore, reduces to computing all the lines that intersect a vertical line segment 
corresponding to e in the dual plane. 

We consider two different approaches to the problem. 



4.2.1 The Range-Searching Approach 

In this approach we attempt to preprocess the set of lines L such that given a vertical 
query segment q , we can report the lines of L that intersect q in 0(log n + w e ) time, 
where w e is the number of intersecting lines. 

A naive approach would be to locate one end-point of the query line-segment in 
the arrangement A of the lines L, and move to the other end-point, recording all 
intersections with the the lines of the arrangement. This turns out to be inefficient 
because the complexity of an individual cell may preclude a constant-time transition 
to an adjacent cell. But this is precisely what we want, and can be achieved at the 
cost of some extra preprocessing. 

In order to reduce the complexity of an individual cell so that we can move from 
one cell to an adjacent one in constant time, we first trapezoidize the cells of the 
arrangement by drawing vertical lines through the vertices of the arrangement. Each 
such line is terminated by the first line in A that it meets, if any at all. This can 
be done using the sweepline technique in 0(n 2 log n) time. The resulting geometric 
structure is shown in Fig. 26. Let us name this planar subdivision A!. 




Figure 2G: New planar subdivision A! 


Though each cell now has the geometric appearance of a trapezium, there can 
be many vertices on its non-parallel sides (see Fig. 27) in which case it would still 
be very complex. Such a trapezium has an upper (resp. lower) anomaly if its upper 
(resp. lower) side consists of more than three vertices. 

For an edge e € A, let A — {ai, <22, .. . ,a u } be the set all vertical edges in A' 
that are incident on e from above and B = {61, 62, ■ • • , b v } be the set of those edges 
that are incident from below. 

For j in [1 . . . v - 1] let {a,-, a i+u . . . , aj} be the subset of A whose end-points on 
e lie strictly between those of bj and bj+i- We extend the edges Oj + i, a i+3 , . . . , a,- + i, 
where l < (j - i ) (see Fig. 28 ) up to an edge of A lying immediately below e. While 
this removes the upper anomalies of all the trapezium that hang on e, it may create 
upper anomalies for the ones that lie immediately below. 

To handle all the regions uniformly, we process A! by levels of the original 
arrangement A from the topmost level to the bottommost one, removing all upper 
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anomalies. And again from the bottom to the top, removing all lower anomalies. 
The data structure underlying the resulting planar subdivision A" is called a hive- 
graph H [Cha86]. 

The time and space complexities of this preprocessing step are in 0(n 2 log n) 
and 0(n 2 ) respectively. 

We preprocess A" for point location in 0(n 2 log n) time [EGS86], We index the 
regions of A 1 ', and store these regions in an array R such that the array index of a 
region is identical to its region index. Along with a region we also store the index 
of the region just above it. Given a query segment q, we first locate its bottom 
end-point in A" in 0(log n) time; we then use its region index to locate it in the 
array R and determine the index of the region just above. We repeat this till we hit 
the top end-point of q. 

We claim the following result. 


Lemma 4.1 A set of n lines can be preprocessed in 0(n 2 log n) time using 0(n 2 ) 
space such that given a vertical query segment we can report all lines that intersect 
it in 0(log n + w e ) query time, where w e is the size of the output. 


It follows from the above lemma that if we are given a set of n lines L and a 
set of 0(n) vertical line segments V, the intersections of the lines in L with the line 
segments in V can be computed in 0(n 2 log n + K) amortized time, where AT is a 
count of all such intersections. 

The complexity of this algorithm is prima facie worse than than that of the 
straight forward algorithm. But we can whittle away a y/n factor from the above 
complexity by preprocessing groups of lines of size yfn. The details are in the lemma 
below. 
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Lemma 4.2 Given a set L of n lines and a set V of 0(n) vertical segments , 
computing the intersections of the vertical segments and the lines in L can be done in 
0(n 3 / 2 log n+K) time using 0{n) space, where K is a count of all such intersections. 


Proof: We divide the problem into m subproblems, each consisting of n/m lines 

and all the vertical segments in V. Solve each subproblem by the algorithm above. 
Summing up the time complexities of the individual subproblems, we find that the 
optimal time complexity of 0(n 3/2 log n + I\) is achieved for m = y/n. The space 
required is 0(n ) since we are considering y/n lines at a time. | 


Knowing the potential witness set of each edge in DT, we can compute the 0 max 
values for each edge e in time proportional to the size of W e . Therefore j3 max value 
for all the edges of DT can be computed in time proportional to K. Thus we have 
the following result. 


Theorem 4.3 The (3 -spectrum of a set of points P for 0 > 1 can be computed in 
0(n 3i/2 log n 4 - K) time using O(n) space, where K is a count of the total number of 
times edges in the DT are witnessed for 0 = oo. 


Similarly, we can obtain the following result. 


Theorem 4.4 The 0-spectrum of a set of n points P for 0 > 0 can be computed in 
0(rc 2 log n) time. 


We note that in this case we have to query 0(n 2 ) vertical line segments and 
hence dividing the lines into groups of size y/h does not help. 
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4.2.2 Sweepline Approach 


We can avoid a complicated data structure like the hive-graph, and at the same time 
take advantage of the fact that we know all query segments beforehand by using a 
sweepline technique. 

In the sweepline status £ we store the lines in L according to a total ordering 
of their intersections. The event queue Q is an union of the set of vertical segments 
V and the intersections of the lines in L, sorted by their x-values. 

The algorithm works as follows. As we sweep the structure, if the next event 
corresponds to an intersection of two lines we update £ by reordering their positions 
in the intersection-order with the sweepline. If the event is a vertical line segment, 
we report all the lines in £ that intersect it. This can be done in 0(log n + w e ) 
time, where w e is the number of lines that intersect the vertical line segment. 

As in the range searching technique we also have the following two theorems 
which follow from the observations above. 

Theorem 4.5 The ft -spectrum of a set of points P for (3 > 1 can be computed in 

0(n 3 / 2 log n + K) time using 0(n) space, where K is a count of the total number of 

times the edges in DT are witnessed for (3 = oo. 

Theorem 4.6 The /3-spectrum of a set of points P for /3 > 0 can be computed in 

0{n 2 log n + K) time using 0(n 2 ) space, where K is a count of the total number of 

times the edges of complete graph are witnessed for f3 = oo. 


4.3 ^-Spectrum 


As we noted in an earlier chapter, there is a natural generalization of the notion of 
a /7-skeleton to that of a A: /7-skeleton obtained by relaxing the emptiness criterion 
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of a 0- lune to let it contain at most k - 1 points. Similarly, the A;0-spectrum is 
a generalization of the 0-spectrum. The kfi - spectrum of a set of n points P is a 
weighted geometric graph G. Each edge e G G assigned with largest 0 value 0 max , 
for which the 0-neighborhood of the edge contain at most k- 1 points. 

To compute the spectrum for 0 > 1, we start with the A;-Delaunay graph (k- 
DG) [SC90], known to be a super graph of a ^0-skeleton. To compute the spectrum 
for 0 values greater than or equal to 0, we start with the complete graph. 

Similar to previous methods we compute the set of witnesses W e for each edge 
e € k-DG. For each witness w e W e compute the 0 value at which w becomes 
witness of the edge e. The j3 max of the edge e is the fcth smallest of these 0 values, 
if \W e \ > & otherwise 0 max = oo. Therefore, we claim the following result. 


Theorem 4.7 The kfi -spectrum of a set P of n points in the plane for 0 > 1, can 
be constructed in 0(kn 3 / 2 log n + K) time using 0(n ) space, where K is a count of 
the total number of times the edges in k-DT are witnessed for 0 = oo. 


Similarly one can establish the following theorem. 


Theorem 4.8 The kfi -spectrum of a set P of n points in the plane for 0 > 0, can 
be constructed in 0(n 2 log n + K) time using 0(n 2 ) space, where K is a count of 
the total number of times the edges of the complete graph are witnessed for 0 = oo. 


4.4 Weighted ^-Spectrum 

Let W{ be the additive weight of the point pi G P. Each weighted point pi can 
be treated as a circle centered at Pi , with radius equal to its weight Wi. The 
weighted distance d w (pi,x ), between a point Pi € P and an arbitrary point x, is 
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d(pi,x) Wi if x is outside and zero otherwise. The weighted distance d w (pi p ) 
between any two points p u pj € P is the distance between B % and B r Let D denote 
the set of circles B{ . 

The additively weighted 0-spectrum is a weighted geometric graph, each edge e 
being weighted with the laigest P value, Pjnaxi for which the /3-neighborhood of e is 
not intersected by any weighted circle. 

In the rest of this section we discuss an algorithm for computing the weighted 
/5-spectrum. 


4.4.1 An Algorithm 

To compute the spectrum for 0 > 1 we start with a super graph which is the 
weighted Delaunay graph (WDG) [Mir93]. For p > 0, we start with the complete 
graph. We discuss only the first case. 

For each edge e € WDG, we compute the set of witnesses for 0 = oo. The 
previous algorithms do not generalize in a straight forward manner, because a circle 
in the primal plane does not have a dual that is a simple object. However, a simple 
observation helps us get around this problem. To explain this, we need the help of 
some notations. The symbol denotes the luneooipuPj)] we will let S k denote 
the smallest infinite strip which encloses the weighted circle B k and is parallel to S e . 

The following observation is very straightforward. 


Observation 4.1 A point p k is witness to an edge e for some value of 0 if and only 
if S k intersects S e . 


The standard point-line dual transformation T maps a non-vertical infinite strip 
into a vertical line segment. Moreover, two parallel non-vertical infinite strips 
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Figure 20: Illustration of Observation 1. 

intersect each other in the primal plane if and only if their corresponding vertical 
line segments overlap in dual plane. In other words, a point p k is witness of an edge 
e for some value of /? if and only if the vertical line segments and T Se overlap 
each other in the dual plane. We use the simpler notations If and I e for T *t and 
Ts e respectively. In terms of this simpler notation the above observation may be 
restated as below: 


Observation 4.2 A point p k is witness to an edge e for some value of /3 if and only 
if I k intersects I*. 


Therefore, to compute the set of witnesses of an edge e € WDG, we find all the 
vertical line segments If that intersect the vertical segment I e in the dual plane. The 
question now is how efficiently can we get all the vertical segments If’ s for each e £ 
WDG ? The following observations provide an insight. Their proofs are immediate 
from the definition of the dual transform. 
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as follows: if <]■> is the top end point of If, delete If from C Eh , since If is below 
Eli- Otheiwise it is the bottom end point of If] insert If in C E j 2 , since EI 2 is in 
If. Similarly update C 1 /.;/,, with respective to the point 93 - There is no change in the 
interval lists, C E h and C E i 3 . 

When the sweepline £ coincides with I c , we report all the intersecting intervals. 
This is done as follows. We find the elementary interval El 6 £ such that the top 
end point of I e is in EL The interval I e intersects all the intervals in C E i- This 
takes care of all intervals, whose top end point is above EL But the top end-points 
of some intersecting intervals may start below that of EI. These can be reported by 
going down from the top end-point to the bottom end-point of I e in the sweepline 
status £. During this, if the top end-point of any interval is encountered we report 
the interval corresponding to this top end-point. 

It is easy to determine the events arising out of the adjacency of a pair of points 
on the sweepline status £. This is the x-value corresponding to the slope of a 
common tangent to the weighted circles (in the primal plane) of the intervals (in 
the dual plane) to which these end-points belong. This can be done in 0(1) time 
for each event (see Fig. 32). 




Figure 32: Finding an event. 
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Lemma 4.9 J he witness set. for all the edges e E WDG can be computed inO(n 2 log n+ 

K) time using ()(n 2 ) space, for fi = oo. 

Proof: At any time, the sweepline status £ contains 2 n end points of n intervals. 

£ can be implemented simply by an array of size 2 n. At x = 0, the sweepline status 
£ can be initialized in 0(n log n) time. Each pair of adjacent end points in £ forms 
an elementary interval. For each elementary interval El € £ we compute a set of 
intervals Chi, and maintain them in a height-balanced tree. All this can be done in 
0(n 2 log n) time. 

The initial events are determined by adjacent end-points in £. We compute 
these events by traversing the sweepline £, from one end to other in 0(n ) time. We 
also add all the vertical segment I,, for all e €WDG. This initializes the event-queue; 
the complexity of this is in ()(n log n). 

Let e be an event which is a change of order of q 2 and q s (see Fig. 31). The cost of 
locating these point in £, and hence the three consecutive elementary intervals Eli , 

EI 2 , and Eh such that \q 2 ,q-f\ = Eh is in 0(log n). Updating Cei 2 can be done in 
0(log n) time, since we need to make at most two deletions and two insertions to 
the list Chi 2 . Inserting the next events of <71 and q%, q-$ and q 2 , and q 2 and <? 4 into 
Q can be done in 0(log n) time. There are at most 0(n 2 ) events since an interval 
can intersect another interval at most constant number of times. Therefore the total 
cost of processing all such events is in 0(n 2 log n) time. 

When the sweepline £ coincides with /«, in 0(log n) time we locate the ele- 
mentary interval El € C such that the top end-point of I e is in EL We report all 
intervals in Cri and travel down the sweepline £ from the top end-point of I e to the 
bottom end-point of J e , reporting intervals whose top end-point are encountered. 
This can be done in 0(w e ) time since there are at most 2w e end points between the 
top and bottom end-points of J e , where w e number of intervals I e intersected. The 
total cost of processing all the vertical segments I e for e €WDG, is in 0{n log n+K) 
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time, where K is a count of the number of times that the edges of the WDG are 
witnessed for (3 = oo. | 

The complexity of the above algorithm is prima facie worse than the brute-force 
algorithm which tests a circle for intersection with each of the 0(n) infinite strips, 
corresponding to each edge of the WDG! For the complexity of this is in O(n), and 
therefore in 0(n 2 ) for n circles. However, if we look at the above algorithm very 
carefully we notice that the problem lies in considering all the intervals in one go. 
The cost of maintaining the sweepline status of n intervals is high and dominates the 
algorithm. So, similar to previous algorithms here also we balance the processing 
time for intervals and finding the witnesses to improve the time complexity. One 
can easily verify the following lemma. 


Lemma 4.10 The set of witnesses for each edge e <E WDG can be computed in 
0{r i 3 / 2 log 7i + K) time using 0(n) space, where K is is a count of the total number 
of times edges in WDG are witnessed for (3 = oo. 

This gives us the following theorem. 


Theorem 4.11 For (5 > 1, the weighted (3-spectrum of P can be constructed in 
0(n 3 / 2 log n + K) time using 0{n) space, where K is is a count of the total number 
of times edges in WDG are witnessed for (3 = oo. 

Similarly, the weighted /3-spectrum of P, for (3 > 0, can be constructed by 
starting with the complete graph on P as the super graph and find set of witnesses 
for each edge. The following theorem is easily proved. 

Theorem 4.12 The weighted (3-spectrum of a set P of n points in the plane, for 
(3 > 0, can be constructed in 0{nf log n + K ) time using 0[n 2 ) space, where K is 
is a count of the total number of times edges are witnessed for (3 = oo. 
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4.5 Conclusion 


In this chapter we have proposed efficient algorithms for computing the /5-spectrum, 
kp . spectrum and, weighted /5-spectrum in the L 2 metric. It would be interesting to 
implement these algorithms and test them on some practical data. 



Chapter 5 


/3-Skeletons in Higher Dimensions 


5.1 Introduction 


In this chapter we extend the definition of a /3-skeleton to higher-dimensions. We 
show that the size of a lune-based /3-skeleton for 0 > 2 in higher dimensions is linear. 
We also present an efficient algorithm for computing a /3-skeleton for /3 > 2. 

The rest of the chapter is organized as follows. In the next section we define the 
two types of /3-skeletons - lune and sphere-based. In the following section, we prove 
an upper bound on the size of a /3-skeleton for /3 > 2. We next present an efficient 
algorithm for computing a /3-skeleton for /? > 2. We conclude in the third and final 
section. 


5.2 Definitions 

A /3-skeleton is a geometric graph, obtained by joining pairs of points whose /?- 
neighborhoods are empty. The neighborhood size depends on the parameter /3, and 
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are of two types. 

The sphere-based ^-neighborhood of a pair of points (p l ,p J ) is defined as follows: 
for /5 > 1, it is the union of all open spheres of radius p\pipj\/2, passing through pi 
and pj. For p € [0, 1], it is the interior of the intersection of all open spheres of 
radius |pipj|/(2/5), passing through p, and pj. 

The lune-based /^-neighborhood of a pair (puPj) of points is defined as follows: 
for (3 > 1, it is the interior of the intersection of two spheres of radius P\p?p]\/2, 
centered at the points (1 - P/2)pi + (P/2)pj and (P/2)pi + (1 - P/2)pj , respectively. 
For P E [0, 1], the lune-based neighborhood is identical with the sphere-based 
neighborhood. 


5.3 Algorithm 

5.3.1 Main Ideas 

The essential technique is to compute a small super graph of a /5-skeleton and then 
prune the edges of this super graph that do not belong to the /5-skeleton. 

To accomplish the first step we use the powerful region approach of Yao [Yao82]. 
This gives us a super graph, called the Geographical Neighborhood Graph (GNG, 
for short). 

In the section below, we discuss an improved version of Yao’s approach. To aid 
the intuition we discuss first the two-dimensional case and then indicate how the 
result can be extended to higher dimensions. 
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5.3.2 The two-dimensional case 

The key idea of this approach is to decompose the neighborhood of a point p t - into 
a finite number of sectors by drawing rays from p*. Each sector is said to be a 
narrow region and is characterized by the fact that if p j and p k are two points in 
a sector then dist(pj,p k ) < max(dist{p i ,p j ),dist{p i ,p k )). This latter inequality can 
be guaranteed by letting the angle of each sector be < 60°. The nice consequence is 
that in each sector only the edges connecting p t - to its nearest neighbors are edges 
of the super graph. 

Our observation is that the upper bound of 60° degrees is a function of the 
parameter value (5 . And though the inequality characterizing a narrow region need 
not hold when we raise the upper bound, the following is true. 



Figure 33: Region approach 


Lemma 5.1 For a given point pi, let the angle of a sector R be 9 = arecos(l/ 

If more than one point in R is nearest to pi, then none of the edges, connecting pi 
to its nearest neighbors is in a (3 -skeleton for f3 > 2. 

Proof: Let the points pj and p k be two nearest neighbors of the point p* in the 

region R (see Fig. 33). Let the angle IpjPiPk = a - The points C { and Cj are the 
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centers of circles whose radius is /3|pipj|/2 and pass through p t and pj respectively, 
as shown in Fig. 34. Elementary geometry shows that the point p k is inside the 
^-neighborhood of p,- and pj\ for it is trivial that CjPj is longer than Cjp k . And the 
fact that pj € /3-neighborhood (p*, p k ) , follows from the upper bound on the angle a. 

I 


Therefore, we divide the region around a point p* by rays that emanate from p*, 
making angles of (m- 1) arccos(l/\/2/3) with the rc-axis, where l<m<6if0<72 
and 1 < m < 5 otherwise. We number the regions between two successive rays 
(m - 1) arccos(l/\/2/3) and m arccos(l/V2/3) by m and label the region bounded 
by these as Rm{Pi)- 

An important consequence of the above lemma is that we get a super graph of 
linear size for /3 > 2. 

To construct the GNG, for each point p, we have to find its nearest neighbors in 
each sector R{. We show how to do this for the region Ri (the remaining regions can 
be handled similarly). We rotate the y-axis by an angle 90° - 9. This is equivalent 
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to transforming each point p into a new point p' by the following transformation. 


Pi' 


1 -cot(6 l ) 


Pi 

. Py' . 


_ 0 l/sin(0) 


.Py . 



Figure 35: Divide the points into cells. 

We sort the new points by their X-coordinates, and divide them into m con- 
secutive subsets of equal size. We sort the points in each subset by their new 
F-coordinates and divide each subset into m consecutive subsets of equal size as 
shown in Fig. 35. For each cell q, we compute the following four attributes: 

X m ax(Ci ) = max{p x '\p € Ci} 

YmaxiCi) = max{Py<\p e Ci} 

Xmin{ci ) = min{p x < \p € q} 

Xmin(ci) = min{p y '\p G c,} 

With respect to p, the m 2 cells are divided into three classes. A cell is in class 
1, if all its points are in Ri{p)- If all its points are outside i?i(p), it is in class 2; 
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otherwise in class 3. To compute the nearest neighbors of p in R^p) we use following 
procedure: cells of class 2 can be ignored. We compute the Voronoi diagram for each 
cell of class 1 in 0(n/rn 2 log(n/m 2 ) time per cell, since the number of points in each 
cell is in 0{n/rn 2 ). We locate p in each cell’s voronoi diagram and find the nearest 
neighbor in each cell in O(log(n/m 2 )) time. We scan all the points in cells of class 
3 and find the nearest neighbor of p that is in R\(p). 

Therefore, the total time required for computing nearest neighborhood in 
for each point is in 0(n log (n/m 2 ) + nm 2 lo g{n/m 2 ) + n 2 /m). Choosing m = 
(n/ logn) 1/3 , the time complexity is seen to be in 0(n 5/3 log 1/3 n ). 


5.3.3 The higher-dimensional case 

To extend the above discussion to higher dimensions, we need the notion of a basis 
and a frame as introduced by Yao [Yao82]. We consider our metric space to be a 
d-dimensional real vector space. If B = bi,b 2 ,...,bd is a basis of this vector space, 
then the convex cone, Conv(B) = {EiLi > 0 for all z} corresponds to what 
we called a sector in the two-dimensional case. A set of bases Bi, B 2 , ..., B k form a 
frame, R, if U Cone(Bi) is the entire vector space. 

The angle of a cone is the maximum angle between any two vectors in this cone 
and the angle of the frame is the maximum angle taken over all the cones. 

Yao proved the interesting result that if ip is any angle in the range 0<ip < ir, 
we can always construct a frame whose angle is less than ip. For our problem, we 
need to construct a frame whose angle is less than arccos(l/ \/2/5). The independence 
of this angle with respect to dimension is a bit of a surprise. 

As in the 2-dimensional case, the following lemma identifies the possible edges 
of a /^-skeleton for /? > 2. 
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Lemma 5.2 For a given point Pi , let the angle of a cone R be bounded by 6 = 
arccos(l / If more than one point in R is nearest to Pi , then none of the edges, 
connecting P i to its nearest neighbors is in a (3- skeleton for (3 > 2. 


Proof: Let the points Pj and Pk be two nearest neighbors of the point Pi in the 

region R. Let h be a plane passing through the points p hPj , <mdp k . The intersection 
of h and the ^-neighborhood^, Pj ) is a 2-dimensional /^-neighborhood. Therefore, 
the point Pk is in /?-neighborhood(p i ,p j ) from Lemma 5.1. Similarly, we' can show 
that the point pj is in the /?-neighborhood(pj,p fc ). | 


As in the 2-dimensional case, we conclude that for fixed d, the size of the /3- 
skeleton is in 0(n) for {3 > 2. 

Given a frame T, the GNG can be computed using the method proposed by 
Yao in [Yao82], Here is a brief sketch of his algorithm. Cnsider a basis Bi. We 
transform the coordinates of our points to this new basis Bi. Sort the points on the 
first coordinate and divide them into m consecutive subsets of equal size. In each 
subset, sort the points on the second coordinate and divide these into m consecutive 
subsets of equal size. Repeat this processes for the remaining d — 2 coordinates. This 
divides the points into m d subset of equal size in 0(dn log n) time. We preprocess 
each subset such that, given a query point p, we can find a nearest neighbor in 
0(log (n/m d )) time. For each point p, its nearest geographical neighbors can be 
computed as follows. Imagine the frame T to be placed at p. Consider a cone R. 
With respect to R , the m d cells are divided into three classes. A cell is in class 1, if 
all its points are in R. If all its points are outside R, it is in class 2; otherwise it is 
in class 3. To compute the nearest neighbors of p in R we use following procedure: ' 
cells of class 2 can be ignored. We locate p in the Voronoi diagram of each cell in 
class 1 and find the nearest neighbor in 0(\og(n/m d )) time for this cell. We scan all 
the points in the cells of class 3 and find a nearest neighbor of p that is in R. This 
process, repeated for all the cones at p and all the input points, gives us the GNG. 
Following the analysis of [Yao82], we can prove the following result. 
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Lemma 5.3 Given a set P of n points, a super graph of the (3-skeleton for (3 > 2 
can be computed in 0(n 2 ~ n{ '' ) (logn) 1_a W) time, where a[d ) = 

5.3.4 Pruning the GNG 

A brute-force method of thinning the super graph to the ^-skeleton is to check the 
status of each edge. This means checking if the lune corresponding to an edge is 
empty or not. Since this requires 0(n ) time and there are 0(n) edges, the complexity 
is in 0(n 2 ). 

To improve on this result, we consider all the spheres that determine the lunes of 
the edges in GNG. The resulting arrangement of spheres divides the space into 0(n d ) 
cells. We can preprocess these hyperspheres, to obtain for each Mace, 0 < k < d, 
the set of lunes that contain this Mace. For each point p, locate the face / in 
the arrangement of hyperspheres and prune the all the lunes that contain /. That 
is, the problem of checking emptiness of lunes is reduced to point location in an 
arrangement of hyperspheres. 

Since an arrangement of hyperplanes is easier to work with, by an inversive 
transform we map the spheres to a set of hyperplanes in d + 1 dimensions. We 
briefly discuss the inversion transformation. 

An inversion is determined with respective to two parameters. These are the 
center and radius of inversion. Consider the planar case, with the center of inversion 
at the origin and radius of inversion l > 0. A point p is mapped into another point p 
such that d(o,p) d{o,p') = 1 In other words, a vector in the direction 6 is mapped to 
another vector in the same direction but with inverted magnitude. Similarly, a circle 
which passes through the center of inversion is mapped to a line which does not pass 
through the center of inversion. Moreover, the interior of the circle is mapped to a 
half-plane as shown in the Fig. 36. 
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Figure 36: Inverse transformation in plane. 


The inversion of our hyperspheres to hyperplanes in a space that is one dimen- 
sional higher is carried out this way. Imagine our spheres to be embedded in a d- 
dimensional hyperplane in d+ 1-space. Let c be a center of inversion, lying outside 
this hyperplane. Consider a hypersphere, S ' , that passes through this center of 
inversion and a given sphere, S , that lies embedded in the d-dimensional hyperplane. 
The d-dimensional hyperplane, that S' inverts to with repect to c corresponds to S. 

Now the problem of point location in the arrangement of hyperspheres in d-space 
is reduced to the problem of point location in the arrangement of hyperplanes in 
(d -f- l)-space. Using the generalized binary search method for higher dimensions, 
proposed by Dobkin and Lipton [DL76], we can preprocess these hyperplanes for 
point location query. We briefly describe the method, first in 2-dimensions and then 
in higher dimensions. 

Given a set of n lines in the plane, these lines are preprocesed for point location. 
The projections of the intersections of these lines on the x-axis are sorted. These 
projected points divide the x-axis into 0(n 2 ) intervals, such that in each interval no 
two given lines intersect. Therefore, for each interval store the relative ordering of 


73 


the given lines. The time required is in 0[n 2 log n + n 2 n log n). The query can be 
answered in O(log n + log n 2 ) time. That is, 0(3 log n) time. 

We proceed similarly in (d+ l)-space. Project the intersections of all pairs of the 
given 7i hyperplanes on to a hyperplane, which is distinct from given hyperplanes. 

In each region in the arrangement of these 0(n 2 ) projected hvperplane, no two 
original hyperplanes intersect. So for each region store the relative ordering of 
these hyperplanes. To locate the region in 0(n 2 ) hyperplanes in d-dimension, we 
recursively solve the problem. 

Let V(n,d + 1) be the time required for preprocessing the n hyperplanes in 
(d + l)-dimension. We get the following recurrence relation for P(n,d): 

V(n, d + 1) = P(n 2 , d.) + 0(n 2d n log n) 

with V(n, 2) = 0(n?n log n). 

Let Q(n, d +1) the time required for answering a query. We have the following 
recurrence relation for Q(n,d): 

Q(n, d + 1) = Q(n 2 , d) -t- 0( log n ) 

with Q(n, 2) = 0(log n). 

The solutions to these rcicurrences allow us to conclude the following. 

Lemma 5.4 A set ofn hyperplanes in (d+1) -dimension can preprocessed in 0(n 2 1 1< 

time so that each subsequent point location query can be answered in 0((2^ +1 ^ 2 + 

(d — 1)) log n) time. 


The complexity of the above algorithm is worse than the brute force algorithm 
for pruning the lunes. However, if we look at the above algorithm carefully we 
notice that the preprocessing time for 0(n) hyperplanes is high and dominates the 
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algorithm, whereas the query time is low. We can improve on the complexity by 
balancing the preprocessing and query times to derive the following result. 

Lemma 5.5 Given a GNG, the nonempty lunes can be eliminated in 
0(n 2 ~ a '^ i 0 g n ) time for (3 > 2, where a'(d ) = ( 2 d+1 — l) -1 ). 

Proof: We divide the problem into 0(n/m ) subproblems, each subproblem con- 

sisting 0(m) hyperplanes and all the input points. For each of these subproblems, we 
preprocess the hyperplanes for point location in 0((n/m) log m) time. The 

total point location time is in 0(n 2 /m log m). Therefore the total time required to 
solve all the 0(n/m) subproblems is in 0((n/m) m 1//a '( d ) log m + n 2 /mlog m). The 
minimum time of 0(n 2 ~ a '^ log n) is achieved for m = n a '^ d \ | 

The only remaining problem is to generate the lune list for each face. This 
can be done with the preprocessing time of hyperplanes [SC91cj. Hence follows the 

theorem. 

Theorem 5.6 A (3-skeleton for (3 > 2 can be computed in (^(^““^(log n) 1 " 0 ^) 

time. 


5.4 Conclusion 


In this chapter, we have generalized the notion of the /5-skeleton to higher dimensions 
and presented an efficient algorithm for computing a lune-based /5-skeleton for (3 > 2 
in any dimension. 


The problem of finding an optimal algorithm for computing ^-skeletons remains 
open. A recent interesting discovery is that in 2-dimension a circle-based /5-skeleton 
is subgraph of a Minimum Weight Triangulation on n points for (3 > 1.7317. It 
would be interesting to know if this result extends to higher dimensions. 



Chapter 6 


Conclusions 


In this thesis, we have attempted to further some of the of the algorithmic aspects 
of the class of geometric graphs known as /3-skeletons. The highlights of our contri- 
bution, in chapterwise order, are as follows. 

In the first chapter, we have have shown how to compute efficiently a /3-skeleton 
for any /3 > 0. The complexity of our algorithm is in 0(n 3 / 2 log n) for /3 > 1 and in 
0(n 5/2 log n) for /3 in in the range [0, 1), for any L p , where 1 < p < oo. In the Li 
and Loo metrics, the algorithms are in 0(n 5/2 log n) for any (3 > 0. 

We have also extended our algorithms to compute generalized /3-skeletons, namely, 
fc/3-skeleton and additively weighted /3-skeleton. The time complexity of computing 
the it/3-skeleton is in 0 (Am 3 / 2 log n + A: 2 n) for /3 > 1 and in 0(n 5 / 2 log n + n 2 k ) for 
/? in the range [0, 1). The time complexity of computing the weighted /3-skeleton is 
in 0(n 3 / 2 log n) for /3 > 1 and in 0(n 5/2 log n) for /3 in the range [0, 1). 

In the metric Too, a /3-lune is simply a rectangle parallel to one of the axes. In 
the second chapter, we have exploited this property to present an optimal output- 
sensitive algorithm for computing a /3-skeleton for /3 > 2. The time required is in 
0(n log n + K), where K is size of the output graph. 
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In the third chapter, we have presented an efficient algorithm for computing 
the /3-spectrum. Our algorithms are in 0(n 3 / 2 log n + K) for p max > 1 and in 
0(n 2 log n + K ) for P max ^ 0, where K is a count of the number of times that the 
edges of the initial (Delaunay and complete respectively) graph are witnessed for 
/3 = oo. 

We have also extended the notion of the /3-spectrum to &/3-spectrum and addi- 
tively weighted /?-spectrum. We have described fast algorithms for computing these 
graphs. The time complexity for computing the ^-spectrum is in 0{kn 3 / 2 log n+K ) 
for /? m ax > 1 and in 0(n 2 log n + K ) time for p max > 0, where K is a count of the 
number of times that the edges of the initial (^-Delaunay and complete respectively) 
graph are witnessed for p — oo. The algorithm for computing the weighted /3- 
spectrum takes 0(n 3 / 2 log n + K) time for p max > 1 and 0(n 2 log n + K ) time 
for /3 max > 0, where if is a count of the number of times that the edges of the 
initial (weighted Delaunay and complete respectively) graph on P are witnessed for 
P = oo. 

We have extended the notion of circle-based and lune-based /3-skeletons to higher 
dimensions. For fixed d, we have shown that the size of the lune based /3-skeleton 
in d dimension is linear for /3 > 2. We have also presented an efficient algorithm for 
computing a lune-based /3-skeleton for /3 > 2. 


6.1 Further Research Problems 

The research reported in this thesis opens up many new questions. 

A rather obvious one is to find more efficient algorithms for computing the 
various geometric graphs that we have discussed in this thesis, viz., /3-skeleton, k0- 
skeleton, weighted /3-skeleton etc. 

In the metrics Li and L computing the /3-skeleton, for p in [1,2), in time lower 
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order than 0{n log n) is an interesting open question. It is also worth investigating 
if the output-sensitive algorithm presented in this thesis for computing the /3-skeleton 
can be generalized to higher dimensions for /3 > 2. 

One really interesting avenue that we haven’t explored is the application of 
/3-skeletons to practical problems. We hope that somebody will take up this applca- 
tion challenge. Like, for example, using /3-skeletons for shape reconstruction and 
identification of real world objects. 

Recently, researchers have dicsovered an intriguing link between circle-based 
/3-skeletons and Minimum Weight Triangulations [Kei94, Yan95, CX96]. For a 
planar point set, a circle-based /3-skeleton is a subgraph of the Minimum Weight 
Triangulation for /3 > 1.17682. It would be worth investigating if this result 
generalises to higher dimensions. 

Another line of work that is tantalizingly open is the problem of recognising 
a /3-skeleton. Some preliminary work in this direction has been done by Bose et 
al [BLL96]. But a vast untapped field lies ahead. 



Appendix A 

Some Examples of /5-Skeletons 




P=1.00 


P=1.25 



p=1.50 
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